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ABSTRACT 

Aims. We test the robustness of published time delays for 1 1 lensed quasars by using two techniques to measure time shifts in their 
light curves. 

Methods. We chose to use two fundamentally different techniques to determine time delays in gravitationally lensed quasars: a 
method based on fitting a numerical model and another one derived from the minimum dispersion method introduced by Pelt and 
collaborators. To analyse our sample in a homogeneous way and avoid bias caused by the choice of the method used, we apply both 
methods to 11 different lensed systems for which delays have been published: JVAS B0218+357, SBS 0909+523, RX J091 1+0551, 
FBQS J0951+2635, HE 1104-1805, PG 1115+080, JVAS B1422+231, SBS 1520+530, CLASS B1600+434, CLASS B1608+656, 
and HE 2149-2745 

Results. Time delays for three double lenses, JVAS B0218+357, HE 1104-1805, and CLASS B1600+434, as well as the quadruply 
lensed quasar CLASS B1608+656 are confirmed within the error bars. We correct the delay for SBS 1520+530. For PG 1115+080 
and RX J091 1+0551, the existence of a second solution on top of the published delay is revealed. The time delays in four systems, 
SBS 0909+523, FBQS J0951+2635, JVAS B1422+231, and HE 2149-2745 prove to be less reliable than previously claimed. 
Conclusions. If we wish to derive an estimate of Hq based on time delays in gravitationally lensed quasars, we need to obtain more 
robust light curves for most of these systems in order to achieve a higher accuracy and robustness on the time delays. 

Key words. Gravitational lensing: strong - Methods: numerical - Galaxies: quasars: individual: JVAS B0218+357, SBS 0909+523, 
RX J0911+0551, FBQS J0951+2635, HE 1104-1805, PG 1115+080, JVAS B1422+231, SBS 1520+530, CLASS B1600+434, 
CLASS B 1608+656, and HE 2149-2745 - Cosmology: cosmological parameters - 



1. Introduction 

The determination of the time delay between different images 
of gravitationally lensed quasars is an important step in differ- 
ent kinds of studies: in deriving Hq, the expansion rate of the 
Universe (e.g.fV uissoz et al. 2008 ), for microlensing studies (e.g. 



Para ficz et al.|2006a , and for detailed studies of the structure of 



Goicoechea|2002||Morgan et al.|2008a|>. 



a lensed quasar (e.g 

However, previous time delay determinations have been far 
from homogeneous, not only because they are based on different 
methods, but also because of their varying levels of reliability. 
Their accuracy depends, among other factors, on the amplitude 
and shape of the quasar's intrinsic variations, the perturbations 
of the light curves by microlensing effects, the photometric error 
bars, the typical time sampling of the monitoring, the total time 
span of the observing campaign, and the intervals between ob- 
serving seasons. Moreover, the published error bars are generally 
internal errors only and the way in which they are determined 
varies from study to study. Unfortunately, once a time delay has 
been published, the value may be used for years without verifi- 
cation, even when the authors of the original article caution the 
reader that the result is not very well-constrained (e.g. Jak obsson] 
|et al.|2005] >. 

Hence, we are of the opinion that it would be useful to re- 
evaluate published time delays in a number of systems using the 
same methods as for all of the delays, as well as for the estimate 
of the error bars. An idea of the robustness of the time delay val- 
ues can thus be obtained, not only internally with our methods, 
but also by comparing our results with published ones. We do 
not claim that our methods are superior to other ones. However, 
a critical reanalysis of published results using two fundamen- 



tally different approaches allows us to sort the results in terms 
of the reliability and independence of the method and to deter- 
mine which lensed systems may be useful for determining Ho. 
Our main purpose is thus to examine whether the published light 
curves allow the determination of reliable time delays. If the an- 
swer is positive, we attempt to estimate realistic error bars, and 
correct some of the published values for small systematic errors. 

For the determination of Ho by means of gravitational lens- 
ing to be competitive with more classical methods (e.g. Riess 
et al.|2009] l, we need to reach at least a comparable accuracy of 
~ 5% in Ho. As time delays from different lensed systems should 
be combined to obtain Ho, we can assume that the error in the 
individual time delays contributes to the statistical error in Hq. 
Since the time delay uncertainties are only one of several sources 
of error in Ho determinations (to be added e.g. to uncertainties in 
the dark matter distribution in the lens), they should in any case 
not exceed 5%. 

Section 2 and 3 present the two methods used for time de- 
lay determination. These methods are applied to the 1 1 lensed 
systems, for which the main results are described in Section 4. 
The lenses of our sample are those for which accurate astrome- 
try has been determined by means of the deconvolution of near- 
infrared Hubble Space Telescope images (Sluse et al, submitted 
to A&A). A summary of these results is presented in Section 5, 
together with our conclusions. 



2. Numerical model fit (NMF) 



We revised and improved the method described in Burud et al. 
( |2001| l. The basic idea can be summarized as follows: for a series 
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of given time delays, the method minimizes the difference be- 
tween the data and a numerical-model light curve with equally 
spaced sampling points, while adjusting the two parameters of 
the difference in magnitude between the light curves and a slope 
that models slow linear microlensing variations. The model is 
smoothed by introducing the convolution of the model curve 
with a gaussian r{t) of full width at half maximum comparable 
to the typical sampling of the observations, and this smoothing 
term is weighted by a Lagrange multiplier A. The function to be 
minimized is: 



x 2 + 



■(r*g)(fdr 
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as used in |Burud et al.| ( |2001) where c4(f,) and o>(f,) are the 
data for image k(k = A, B) with the associated error bar, g(tj) the 
model curve, At the time delay, and Am and a the parameters 
representing the difference in magnitude and slope between the 
light curves. 

The optimal time delay is the one that minimizes the reduced 
X 2 red between the model and the data points. It is important to 
insist on the difference between^ 2 as defined in Eq.[2]and on the 



other hand the reduced 
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in which x 2 is divided by the number of data points in com- 
morQbetween the light curves of the quasar images for a given 
time delay. Indeed, the longer the time delay one tests, the fewer 
points these light curves have in common, which tends to reduce 
the x 2 and result in a bias towards longer time delays, hence the 
use of x 2 ed to avoid this bias. 

A second important difference from the original version of 
the method is technical: for computational reasons, the length of 
the model curve should be a power of two, which in some cases 
proves to be too long in comparison to the data, thus falsifying 
the balance between data and smoothing terms. In the original 
version the smoothing term was applied to the full length of the 
model. We adapted the program in such a way that the part of the 
model that is only needed to complete the length until the next 
power of two, is no longer taken into account in the minimization 
process. In this way, the method becomes independent of the 
number of data points in the light curve, which was not the case 
in its original form. 

This method has not only been implemented for two light 
curves, but also for deriving time delays from three and four 
light curves simultaneously. The strength of this simultaneous 
approach lies not only in the improved constraints on the model, 
but also in that we assume the coherence between pairs of time 



1 By this, we mean the data points lying in the time span for which 
data for all the lensed images are available after shifting the light curve 
for the assumed time delay. The number of data points in this common 
time span is not necessarily the same for the light curve of every lensed 
image, hence the use of Na and N B . 



delays, differences in magnitudes, and the slope parameter val- 
ues. 

The robustness of the measured time delay is tested in two 
ways. First of all, we iteratively attempt to find the three pa- 
rameters of the model light curve: the spacing of the model 
curve's sampling points, the range of the smoothing term, and 
the Lagrange multiplier. The results should be independent of 
these parameters as long as we remain in a certain range adapted 
to the data. 

In a second step, we wish to test the influence of each indi- 
vidual point of the light curve on the time delay. This is achieved 
by means of a classical jackknife test: for a light curve consist- 
ing of N data points, we recalculate N times the time delay in 
the light curve of N-l data by successively leaving out one data 
point at a time. Time delays should not change drastically be- 
cause of the removal of a single point from the light curve. If 
they do, we know which data point is responsible for the change 
and we can have a closer look at it. 

Errors are calculated by means of Monte Carlo simulations. 
Normally distributed random errors with the appropriate stan- 
dard deviation are added to the model light curve and the time 
delay is redetermined. We note that errors are not added to the 
data as they already contain the observed error, so adding another 
would bias the results. The model, to which the measurement er- 
rors are added, is assumed to provide a more accurate description 
of the real light curve of the quasar than the data. This procedure 
is repeated at least 1000 times, preferably on different combina- 
tions of smoothing parameters. The mean value of the time delay 
distribution that we obtain is considered to be the final time de- 
lay and its dispersion represents the 1 <x error bar. When we have 
a markedly asymmetrical distribution, we take its mode as the fi- 
nal time delay and use the 68% confidence intervals to obtain 
error bars. In this paper, all quoted uncertainties are 1 cr error 
bars except where mentioned explicitly. 

The advantages of this method are manifold. First, none of 
the light curves is taken as a reference curve; they are all treated 
on an equal basis. Second, a model light curve is obtained for 
the intrinsic variations in the quasar, which is also the case for 
the polynomial fit method described by Kochanek et al. ( 2006|l, 
but not for the minimum dispersion method developed by Pelt 
et al. ( 19961. This is important when calculating the error bars, 
thus avoiding adding random errors to the data. Finally, since 
the model is purely numerical, no assumption is made about the 
quasar's intrinsic light curve, except that it is sufficiently smooth, 
and we only interpolate the model, never the data. 



3. Minimum dispersion (MD) 

The second method we use is derived from the minimum disper- 
sion method b y|Pelt et al.|(|1996|l w ith a number of adjustments 
as describe d in|Courbin et al.| ( |2010| l. The main improvements to 
the original |Pe It et al.| ( |1996| l method consist in: 

1 . No light curve is taken as a reference, they are all treated on 
an equal basis; 

2. A flexible modelling of microlensing by polynomials up to 
third order, per light curve or per observing season. 

Since no model light curve is constructed, computation time 
is a lot shorter than for the NMF method. By using two methods 
based on completely different principles, we are able to check 
whether the derived time delays are independent of the method, 
thus testing their robustness (i.e. independence of the particular 
way in which the data are analysed). 
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4. Application to 11 lensed quasars 

We now present the main results of our time delay analysis for 
each of the published light curves of 11 gravitationally lensed 
quasars. 

- JVAS B0218+357 We have used the data set published by 
Cohen et al.| ( |2000| >, consisting of 5 1 flux density measure- 
ments at 8.4 GHz and 15 GHz. The authors obtained a time 
delay At A s = 10. days, where A is the leading im- 
age, thus confirming independently two values published 
earlier by Big gs et al.| ([1999 ) and |Corbett et at] (1996) of 
AtAB = 10.5 ± 0.4 days and At AB = 12 + 3 days, respectively. 
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Fig. 1: Light curves at 8.4 GHz and 15 GHz for JVAS 
B0218+357 after transforming the flux density measurements 
into magnitudes. The B curve has been shifted by 1 magnitude 
for clarity. 



After transforming the flux densities onto a logarithmic scale 



as shown in Fig. [T 
GHz and 15 GHz 



we applied the NMF method to the 8.4 
ight curves. Using the entire light curve 
did not give a clear and unique solution. The jackknife test 
shows that certain data points can change the value of the 
time delay. After eliminating three of these points, the 9th, 
12th and 35th, from the 8.4GHz light curve, and choosing 
appropriate smoothing parameters, we obtain a time delay 
AtAB = 9.8+q g days at 68% confidence level. The larger er- 
ror bars for higher values of the time delay are due to a sec- 



ondary peak in the histogram (see Fig. [2} around AtAB ~ 14 
days. 




Time Delay (days] 



Fig. 2: Histogram of 1000 runs of the NMF method for the 8.4 
GHz data light curve of JVAS B02 18+357 leaving out three de- 
viating points. One day gaps in the histogram are artefacts due 
to the quasi-periodicity of the data. 



Taking into account all points of the 15 GHz light curve 
provided a comparable value of the time delay of At A s = 
1 1 . l^l j, even though we noted that the importance of the 
secondary peak around AtAB ~ 14 days was significantly 
lower after we had eliminated outlying points in both the A 
and B curve, in the same way as for the 8.4 GHz curve. This 
suggests that the secondary peak around AtAB ~ 14 days is 
probably caused by artefacts in the data, hence we can con- 
firm with confidence the previously published results: com- 
bining the values based on the 8.4 GHz and 15 GHz light 
curves gives a time delay of At A B = 9.9*q° days. 
The MD method confirms the strong influence of these devi- 
ating points: the secondary peak around AtA B ~ 14 days even 
completely disappears when they are removed from both the 
8.4 GHz and the 15 GHz light curves. The 8.4 GHz curve 
gives a time delay of At AB = 12.6 ±2.9 days, and the 15 GHz 
data lead to At A s = 1 1.0 ± 3.5 days, which gives a combined 
result of AtA B = 11.8 ± 2.3 days, all in agreement with the 
above-mentioned values. 

Even if all of these time delay values for this object are in 
agreement with each other, the data do not allow a precision 
of the order of 5% in the delays, which would be necessary 
for a useful estimate of Hq. 

SBS 0909+523 We us ed the data set published by 
p008l , 



Goicoechea et al. 



which contains 78 data points 
spread over two observing seasons. Their analysis leads to 
a time delay AtsA =49 + 6 days where B is the leading im- 
age, confirming the previously reported delay AtsA - 45+j 1 



of |Ullan et aL] l|2006). 
The NMF method, when applied to the entire light curve, 
gives a delay AtsA ~ 47 days, as displayed in Fig. [3] which 
is within the error bars of the previously published delay. On 
closer inspection however, we note that this delay strongly 
depends on two points that are outside the general trend of 
the lightcurve for image B and fall right at the end of the 
time interval covered by the A data points for this time delay 
value: the 63rd and the 64th data points. Recalculating the 
delay while omitting these two points gives a different result 
of AtBA ~ 40 days or even lower values, which is not within 
the published ranges. The same happens if we only take into 
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account the second observing season, which is the longer 
one: the delay then shortens to AtgA ~ 40 days. The param- 
eter modelling slow linear microlensing is also significantly 
smaller in this case. Visually, both results, with and with- 
out the two problematic points, are acceptable. Nevertheless, 
although both values have proven to be independent of the 
two smoothing parameters, the NMF method is sensitive to 
the addition of normally distributed random errors with the 
appropriate standard deviation at each point of the model 
curve. This is because the dispersion in the data points is too 
small compared to the published error bars. That we obtain 
X 2 red <s 1 also highlights some possible problems in the data 
reduction or analysis. 



3600 3650 3700 
JD- 2450000 (days) 



Fig. 3: Light curves of SBS0909+523: A is shifted by a time de- 
lay A, = 47 days and a difference in magnitude Am = -0.6656. 
The slope parameter a = 6.0587 • 10~ 5 corresponds to a model of 
slow linear microlensing. The encircled points are the 63rd and 
64th observations that were omitted from later tests. 



The MD method gives similar results. When using all data 
points, two possible delays can be seen, depending on the 
way microlensing is modelled: AtBA ~ 49 and AtgA ~ 36. 
When the two aforementioned data points are left out, we 
only find At BA ~ 36, independently of how microlensing is 
handled. 

In all cases, with or without these points, leaving more or 
less freedom for the microlensing parameters, the large pho- 
tometric error bars result in very large error bars in the time 
delay when adding normally distributed random errors to 
the light curves, so that delays ranging from At BA ~ 27 to 
AtBA ~ 71 are not excluded at a 1 cr level. 
We conclude that this light curve does not allow a reliable 
determination of the time delay. To determine whether these 
two data points that do not follow the general trend are due to 
genuine quasar variations and thus crucial for the time delay 
determinations or whether in contrast, they are affected by 
large errors and contaminate the published results, we will 
need new observations and an independent light curve. 
- RX J0911+0551 Data for this quadruply lensed quasar were 
made available by |Paraficz et al.|([20 06b ), but had been previ- 
ously treated and analysed by[Burud ( 200T]l and |Hjorfh et al.| 
( 2002) 1, who proposed time delays of At BA = 150 + 6 days 
and AtsA = 146 + 4 days respectively, where B is the leading 
image of the system and A the sum of the close components 
Al, A2, and A3. 

Using all data except the first point, which has too strong 
an influence on our slope parameter because of its isolation 



50400 50600 50800 51000 51200 51400 51600 51800 52000 52200 
JD - 2400000.5 (days| 

Fig. 4: Light curves of RXJ091 1+0551: A, which is the sum of 
the close components Al, A2 and A3, has been shifted by one 
magnitude for clarity. 



as can be seen in Fig. |4] we can at first sight confirm the 
published delays: the NMF method gives At B A = 150 + 2.6 
days and the MD method results in At B A = 147.4 + 4.6 
days. However, the histogram in Fig. [5] shows a secondary 
peak at AtpA ~ 157 days. Investigating this peak further, 
we come to the conclusion that some points have a very 
strong influence on the delay: the first observing season, 
and especially the first ten points of the light curve, indi- 
cate a shorter time delay. According to Burud (2001 ), these 
points were added to supplement the regular monitoring data 
of the Nordic Optical Telescope. However, the first three 
points of the regular NOT monitoring in the B curve are simi- 
larly crucial. Omitting these three points leads to larger error 
bars of At B A = 151.6 + 7.0 days using the NMF method. 
Finally, recalculating the time delay in the regular moni- 
toring data only and without the first three points in the B 
curve, gives AtsA = 159 ± 2.4 days with the NMF method. 
The MD method results in this case in a histogram with 
two gaussian peaks, one around AtsA ~ 146 days and one 
around AtgA ~ 157 days, implying a mean time delay of 
AtBA = 151.4 ± 6.7 days. Only a new and independent light 
curve of similar length could tell us with more confidence 
which of these values is correct and which is possibly biased 
(e.g. by microlensing). 





Time Delay (days] 



Fig. 5: Sum of three histograms of 1000 runs each for 
RXJ09 11+05, using three different combinations of smoothing 
parameters. 
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- FBQS J095 1+263 5 We used the data s et containing 58 
points published by Paraficz et al. ( 2006b i and presented in 
Fig.0 
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Fig. 6: Light curves of FBQ095 1+2635. The B curve has been 
shifted by 0.8 magnitude for clarity. 



Jakobsson et al. (2005 1 published a time delay At AB = 16 + 2 



days, a result that is only based on the last 38 points of the 
light curve when the system had been observed more inten- 
sively. They found various possible time delays according to 
the method, the smoothing, and the data points included, so 
we performed the same tests. 

We can confirm that the time delay is very sensitive to the 
choice of smoothing parameters in the NMF method, es- 
pecially when using the entire light curve, but is still more 
sensitive to the data points used: leaving out a single point 
completely changes the time delay. We calculated time de- 
lays in the light curve using between 55 and 58 data points 
and we found delays ranging from AtAB = 14.2 ± 4.5 days 
to AtAB = 26.3 + 4.7 days. Taking into account three pos- 
sible smoothing combinations and four sets of data (leaving 
out one more data point in each set) leads to a combined 
histogram of 12000 Monte Carlo simulations, as shown in 
Fig. [7] It is clear that a mean value with error bars At AB = 
20.1 + 7.2 days is not of any scientific use: the error bars are 
too large relative to the time delay. Moreover, the histogram 
is quite different from a normal distribution. There is no sig- 
nificant concentration of the results, which would allow the 
determination of a meaningful mode, independently of the 
chosen binning. One can see that different time delays are 
possible and can be divided in two groups: shorter values of 
~ 10.5, ~ 15, and ~ 18.5 days, and longer values of ~ 26.5 
and 30.5 days. 

When using only the third observing season, which is more 
finely sampled, a relatively stable time delay AtAB = 18.8 ± 
4.5 is measured, but once a single point is left out (for exam- 
ple the 19th point of this third season, a point that deviates 
from the general trend in spite of a small error bar), the result 
completely changes towards longer values (AtAB = 25.0+4.9 
days) and becomes sensitive to smoothing. As the measured 
time delay should not depend on the presence or absence of 
a single point, we can only conclude that this light curve, 
even if it consists of three observing seasons, does not allow 
a precise determination of this delay. 

The MD method entirely confirms the large uncertainty in 
this time delay: using all data points we find a time delay of 
AtAB = 21.5 + 6.8 days, whereas the third season only leads 
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Fig. 7: Sum of twelve histograms of 1000 runs each for 
FBQ095 1+2635, using three different combinations of smooth- 
ing parameters for four sets of data consisting of 58, 57, 56, and 
55 data points. 



to AtA B = 19.6+7.6 days, with peaks in the histogram around 
At AB ~ 12 and At A B ~ 28 days. 

According to Schechter et al.| ( |1998[ > and Jakobsson et al. 
(2005), there are spectroscopic indications of possible mi- 
crolensing, so this might explain the difficulty in constrain- 
ing the time delay for this system. Longer and more finely 
sampled light curves might help us to disentangle both ef- 
fects. However, at the present stage, we can conclude that 
this system is probably not suitable for a time delay analysis 
HE 1104-1805 We used the data published by ~ 



Poindexter 

[eTaTI p00"7l >, which combine th eir own SMARTS R-b and 
data with Wise R-band da ta from|Ofek & Maoz] (|2003] l and 
OGLE V-band data from |Wyrzykowski et al.| ( |206"3"] k The 
three data sets are shown in Fig. |8J Table [1J lists the four 
time-delay values published for HE1 104-1805, where B is 
the leading image. 
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Fig. 8: Light curves of HE1 104-1805, combining the OGLE V- 
band data, the Wise R-band data and the SMARTS R-band data. 



We performed tests with both methods on different combi- 
nations of the data using all telescopes or only one or two 
of them. Unfortunately, the results seem to be sensitive to 
this choice, as they are to the way in which microlensing is 
treated: both the OGLE and Wise data sets were analysed 
to find a time delay Ats A ~ 157 days, whereas SMARTS 
data converge to a higher value of At BA ~ 161 days or more 
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[Morgan et"aT]j20U«a[ ) 



Table 1: Published time delays for HE1 104-1805. 



Poindexter eFaLlp007) 's 



as is shown in Fig. [9] In addition. 

smaller value is recovered with the MD method when in- 
cluding OGLE and Wise data but only for some ways of 
modelling microlensing. We therefore conclude that we can 
neither make a decisive choice between the published val- 
ues, nor improve their error bars, which are large enough to 
overlap. 
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Fig. 10: Light curves of PG 1 1 15+080, where the A curve is the 
sum of the Al and A2 component. The A and the B curve have 
both been shifted by 1.9 and -0.3 magnitude, respectively. 
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Fig. 9: Sum of three histograms for HE1 104-1805, using three 
different combinations of smoothing parameters combining only 
OGLE and SMARTS data. OGLE data point to a time delay 
Mba ~ 157 days, whereas SMARTS data converge to a longer 
value of At BA ~ 161 days. 



- PG 11 15+080 We used the data taken by |Schechter et al.| 
( [1997) . They published a time delay At CB = 23.7 + 3.4 days 



between the leading curve C and curve B and estimated the 
delay between C and the sum of Al and A2 at AtcA ~ 9.4 
days. |Barkana| ( [1997} used the same data, which are shown 
in Fig. 10 but a different method to determine the delays. His 
value of AtcB = 25.0^'g is compatible with Schechter's one. 
Morgan et al. ( |2008b| l published new optical light curves for 



this quadruply imaged quasar in order to study microlens- 
ing in the system. Unfortunately, these light curves cannot 
be used to determine a time delay independently because of 
the clear lack of features in the variability of the quasar and 
the inconsistency of the individual error bars relative to the 
dispersion in the data. 

A first series of tests on Schechter's data with the NMF 
method led to time delays of AtcA ~ 15 days and AtcB ~ 
20.8 days with a minor secondary peak around At CB ~ 
23.8. We then corrected the published data for the existing 
photometric correlation between the quasar images and the 
two stars used as photometric references, as mentioned by 
Barkana ( 1997] >. This caused the shorter time delay to shift 
either towards Ate a ~ 1 1 days or AtcA ~ 16 days, and trans- 
formed the longer delay into two nearly equally possible re- 
sults of Atc B ~ 20.8 days or AtcB ~ 23.8, as indicated by the 
two main peaks in the histogram in Fig.[TT] Adding observa- 



tional errors to the model light curves and taking into account 
four different ways of smoothing results in Atr A = 1 1 .7 ± 2.2 
(see Fig. [H]) and At CB = 23.8^® (see Fig.pt 




Fig. 11: Sum of four histograms of 1000 runs each for AtcB in 
PG1 115+080, using four different combinations of smoothing 
parameters. 




Time Delay (days] 



Fig. 12: Sum of four histograms of 1000 runs each for AtcA in 
PG1 115+080, using four different combinations of smoothing 
parameters. 
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The MD method confirms a delay of AtcB ~ 20.0 days, but 
finds a second solution around AtcB ~ 12.0 days, which re- 
sults in AtcB - 17.9 + 6.9. The value for AtcA is also shorter 
than the one obtained with the other methods namely AtcA = 
7.6 ± 3.9 days and has larger error bars. Unfortunately, the 
length and quality of this light curve do not allow one to 
choose between the possible time delays that differ accord- 
ing to the method used, but are generally lower than pub- 
lished values. 

- JVAS B1422+231 For this quadruply lensed quasar, we used 
the data published by Patnaik & Narasimha (2001), consist- 
ing of flux density measurements at two frequencies, 8.4 and 
15 GHz. Their results for the time delays were based only on 
the 15 GHz data without image D, which is too faint. These 
data are shown in Fig. [13] They obtained AtsA = 1.5 ± 1.4 
days, AtAc = 7.6 + 2.5, and Atgc - 8.2 ± 2.0 days when 
comparing the curves in pairs. 
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Fig. 14: Sum of four histograms for At B A of 1000 runs each for 
JVAS B 1422+231, using four different combinations of smooth- 
ing parameters. 
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Fig. 13: Light curves of JVAS B 1422+231 for the A, B, and C 
components after transforming the flux density measurements at 
the 15 GHz frequency into magnitudes. The C curve has been 
shifted by half a magnitude. 



The NMF method allows time delays to be tested for the 
three light curves simultaneously, thus imposes coherence 
on the results. Given the error bars in the published results, 
which are large compared to the time delays, we can confirm 
the results here, but we emphasize that they include two dis- 
tinct groups of solutions between which we cannot decide 
based on the actual light curves: for the shortest delay, we 
either have At BA ~ 1.0 day or At BA ~ 2.0 days, as shown 
in Fig. [14] The choice between both solutions is sensitive to 
the smoothing parameters: the importance of the first group 
AtgA ~ 1-0 day is lower, and even disappears completely, 
with greater smoothing. 

For At BC , the situation is similar but even less clear: low 
smoothing parameters lead to a range of possible solutions 
between At B c ~ 6 and At B c ~ 10 days, which are all within 
the error bars of the published results. However, Monte 
Carlo simulations of reconstructed light curves with higher 
smoothing parameters give a time delay At B c = 10.8 +1.5 
days with a secondary peak around At BC ~ 8 days. 
The MD method gives a completely different result: it con- 
verges towards time delays that invert the BAC-order into 
CAB but with error bars large enough not to exclude the 
BAC order of At BA = -1.6 + 2.1 days, At AC = -0.8 ± 2.9, 
and At BC = -2.4 + 2.7 days. 



New observations are clearly necessary to reduce the uncer- 
tainties in both the different solutions and the error bars that 
are too large in comparison with the time delays to be useful 
to any further analysis. 

SBS 1520+530 Two data sets exist for this doubly lensed 
quasar: the set made available by Burud et al. (|2002b[) and 
the one published by Gaynul lina et al.| ( |2005a[ ). |Burud et al.| 
(2002c) were the first to publish a time delay for this sys- 
tem AtAB - 128 ± 3 days, where A is the leading image, 
or AtAB - 130+3 days when using the iterative version of 
the method ( |Burud et aL][200T) . |Gaynullina et al.| ( |2005b| l 
used an independent data set and found four possible time 
delays, of which the one with the largest statistical weight, 
At AB = 130.5 + 2.9, is perfectly consistent with the previ- 
ously published time delays. 

Even if the light curve based on |Gaynullina et al.| ( |2005a 
data contains more than twice as many data points as Burud 



et al. ( 2002b[ )'s older light curve, we decided not to use it 



because of the lack of overlapping data between the A and 
the B curves of the quasar after shifting the B curve for the 



time delay, as can be seen in Fig. 15 
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Fig. 15: Ligh t curves for SBS 1520+530 based on |Gaynullina 
et al.| ( |2005"a] rs data after shifting the B curve by At AB = 130.5 
days. Hardly any data points overlap between the light curves 
for images A and B. 
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Applying different tests to |Burud et al.| ( |2002c| l's light curves, 
shown in Fig. 16 using the NMF method led to a time delay 
for which the error bars overlap with the published value: 
AtAB - 126.9 ± 2.3. That the delay is slightly shorter than 



Burud et al. ( 2002c I's value can be explained by our use 
of the reduced X~Td instead of the x 2 , the latter implying 
that longer delays are the more likely ones, as explained 
in Section [2] This effect was also noted using the iterative 
method. 
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Fig. 16: Light curves for SBS 1520+530 based on |Burud et al. 
d2002cl)'s data. 




Fig. 17: Sum of three histograms of 1000 runs each for 
SBS 1520+530, using three different combinations of smoothing 
parameters. 



The MD method yields a comparable time delay of AtAB = 
124.6 ±3.6 days and confirms the shape of the histogram: the 
highest peak value at At AB ~ 125 days and a clear secondary 
peak at AtAB ~ 127.5 days. Combining the values from both 
methods implies that At A B = 125.8 + 2.1 days. 
CLASS B1600+434 Data, as shown in Fig. ~ 



18 were made 



available by Paraficz et al. (2006b) but had been treated 
and analysed by |Burud et al.| ( |2"000| ), who published a final 
time delay At AB — 5 1 ± 4 days, where A is the leading im- 
age, consistent with the time delay At A B = 47^ 2 days from 
Koopmans et al. (2000 ) based on radio data. 
The NMF method leads to a time delay At AB = 46.6 ±1.1 
days, but the histogram in Fig. 19 clearly shows that we 
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Fig. 18: Light curves for CLASS B1600+434: 41 data points 
spread over nearly two years. 



two distinct values of ~ 46 or ~ 48 days, so we prefer to 
speak of a delay of either At A B = 45.6+H days (68% error) 



or At AB = 45.6 + _ 2 Q l days (95% error). 




Time Delay (days) 



Fig. 19: Sum of seven histograms of 1000 runs each for 
B 1600+434, using seven different combinations of smoothing 
parameters. 



These results are in marginal disagreement with the final de- 
lay pro posed by |Burud et al.| ( [2000| . However, |Burud et"aL] 
(2000)'s final result is an average of four time delays, each of 
them calculated with a different method. Two of these meth- 
ods also inferred a value of around ~ 48 days. 
We could identify at least three explanations of our lower 
value and smaller error bars in comparison with Buru d et al 
(2000)'s time delay. The first one is the same again as for 
SBS 1520+530: our use of the reduced x 



2 red (see formula 



cannot use the mean as the final value. The histogram has 



[3] in Section [2]) instead of the x z , the latter introducing a 
bias towards longer delays. The second reason is the tech- 
nical issue concerning the length of the model curve as ex- 
plained in Section 2, which was found to be crucial for the 
time delay of this system. Finally, we observed that higher 
values of the Lagrange multiplier weighting the smoothing 
term seemed to lead to longer time delay values, which dis- 
appeared with lower smoothing. Taking into account these 
three adjustments, nearly all values around ~ 51 days disap- 
pear from the histogram. 

This is not the case for the MD method, which explains the 
slightly longer value of the time delay: At A B = 49.0 ±1.2 



8 



E.Eulaers and P. Magain: Time delays for 1 1 gravitationally lensed quasars revisited 



days. Combining these results gives a delay of AtA B = 47.8 ± 
1 .2 days. Even if these error bars imply that the time delay is 
very tightly constrained, we emphasize that the delay mea- 
surement is only based on 41 data points spread over nearly 
two years, which gives a relatively high weight to every sin- 
gle data point. When adding random errors, neither of the 
two methods leads to a histogram with a gaussian shape. A 
more finely sampled light curve might remedy this situation. 
- CLASS B1608+656 

Light curves for this quadruply lensed system were first anal- 
ysed by Fassnacht et al. (1999) and subsequently improved 
in Fassnac ht et al.| ( |2002| l by adding more data. Using three 
observing seasons, they published time delays of At B A = 
31.5% At BC = 36.0 ± 1.5, and At BD = 77.0+2 days. Their 
analysis is based on a simultaneous fit to data from the three 
seasons but treats the curves only in pairs. 




Time Delay (days) 



Fig. 21: Histograms for the three time delays in B1608+656, us- 
ing two different combinations of smoothing parameters for all 
seasons and three out of the four curves simultaneously. 
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Fig. 20: Light curves for CLASS B1608+656. The B curve has 
been shifted by half a magnitude for clarity. 



We performed several tests on these light curves, which are 
shown in Fig. 20 first taking into account only the first and 
the third season separately (the second season not present- 
ing useful structure), then all data for the three seasons si- 
multaneously, using the three and four curve version of our 
method as described in Section[2] This enables us to impose 
coherence between the pairs of time delays, which was not 
done by |Fassnacht et al.] ( |2002| l. The results, as illustrated in 
Fig. 21 confirm the previously published values, within the 
error bars, of At BA = 30.2 ± 0.9, At BC = 36.2 + 1.1, and 
At BD = 76.9 + 2.3 days. 

One particularity deserves more attention: the value for At B A 
changes slightly according to the seasons and the curves 
considered simultaneously. When leaving out the second 
season data (featureless), At B A systematically converges to- 
wards At B A - 33.5 + 1.5 days, as shown in Fig. [22] This 
is consistent with Fassnacht et al. (2002]), who already men- 
tioned a time delay of AtAc ~ 2.5 days. Even if this slight 
difference is probably due to microlensing and should be in- 
vestigated in more detail, we chose to retain the final value, 
which is the one based on the use of all data. 
The MD method entirely confirms these results of At B A = 
32.9 + 2.9, At BC = 35.2 + 2.5, and At BD = 78.0 + 3.7 days 
with another indication for Af^c ~ 2.5 days. Combining both 
methods results in the time delays of At B A = 31.6 ± 1.5, 
At BC = 35.7 ± 1.4, and At BD = 77.5 ± 2.2 days. 
- HE 2149-2745 We reanalysed the data set made available by 
|Burud et al.| ( [20 02b). These data consist of two light curves, 



Fig. 22: Sum of four histograms of 1000 runs each for At B A in 
B 1608+656, using two different combinations of smoothing pa- 
rameters on the first and the third season. 



one in the V-band and one in the /-band, as shown in Fig. [23] 
|Burud et al[] ( |2002a| l published a time delay At AB = 103 ± 12 
days where A is the leading image. This delay is based on the 
V-band data, but, according to |Burud et al. ( 2002a| l, agrees 
with the /-band data. 

Our tests, based on the light curves as such, and using both 
methods, clearly reveal two possible delays: one around ~ 
70-85 days and one around ~ 100-1 10 days. Unfortunately, 
the light curves of images A and B show little structure and 
hardly overlap, except for some points in the second sea- 
son, when shifting them for a delay of more than 100 days, 
especially in the /-band, which makes it very difficult to 
choose between the two possibilities. Moreover, once we add 
random errors to the model light curve and perform Monte 
Carlo simulations, we only obtain a forest of small peaks, 
spread over the entire tested range of 50 - 140 days, instead 
of a gaussian distribution around one or two central peaks, 
demonstrating that these results are highly unstable. Leaving 
out two outlying data points in the B curve only slightly im- 
proves the situation: within the forest of peaks in Fig. [24] 
those in the range 75 - 85 seem to be slightly more impor- 
tant than those over 100 days. Nevertheless, we cannot derive 
a reliable time delay from these data sets for this system. 
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Fig. 23: Light curves of HE 2149-2745 in the V-band and the I- 
band. The B curve has been shifted by one magnitude for clarity. 
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Fig. 24: Sum of six histograms of 1000 iterations each for HE 
2149-2745 leaving out two data points, using six different com- 
binations of smoothing parameters. 



5. Conclusions 

We have presented an improved numerical method, the numer- 
ical model fit method, to calculate time delays in doubly or 
quadmply lensed quasar systems, and applied it to 11 systems 
for which time delays had been published previously. This al- 
lowed the validity of these time delay values to be evaluated in a 
coherent way. The use of a minimum dispersion method allowed 



us to check the independence of the results from the method. The 
results are summarized in Table [2] 

We caution that some published time delay values should 
be interpreted with care: even if we have been able to confirm 
some values (time delays for JVAS B0218+357, HE 1 104-1805, 
CLASS B 1600+434, and for the three delays in the quadruply 
lensed quasar CLASS B1608+656) and give an improved value 
for one system, SBS 1520+530, many of the published time de- 
lays considered in our analysis have proven to be be unreliable 
for various reasons: the analysis is either too dependent on some 
data points, leads to multiple solutions, is sensitive to the addi- 
tion of random errors, or is incoherent between the two different 
methods used. 

Given the accuracy that is needed for time delays to be use- 
ful to further studies, we note that it will be necessary to perform 
long-term monitoring programs on dedicated telescopes to ob- 
tain high-quality light curves of lensed quasars, not only for new 
systems but also for the majority of the lenses in this sample, 
for which the time delay has been considered to be known. The 
COSMOGRAIL collaboration, which has been observing over 
20 lensed systems for several years now, will soon be improving 
the time delay values for some of these systems for which the 
accuracy is unsatisfactory. 
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System 



Our Result / Comments 



Published Time Delay (days) Reference 

Cohen et al.|{2000} 
Corbettetal\|(|lWo| 
Biggs e t al. ( 1999 ) 
Goicoechea et al. ( 2008 ) 
Lilian et al. ( 2006 ) 



JVAS B02 18+357 

SBS 0909+523 
RXJ091 1+0551 



At A 



. Q Q+4.U 
y - y -0.9 



At AB = 11.8 + 2.3 
unreliable 

2 solutions: 
At BA ~ 146 or ~ 157 



AtAB 

At AB = 12 ±3 
At AB = 10-5 ± 0.4 
At BA =49 + 6 
At BA =45 + _\ l 
At BA = 150 + 6 
AtgA = 146 + 4 



FBQS J095 1+2635 



unreliable 



At A 



16 + 2 



HE 1104-1805 



impossible to distinguish 
but identical within error bars 



AtgA = 152 
At BA 
At BA = 157 
At BA 



+2.8 
3.0 

161+7 



10 



PG 1115+080 



dependent on method 



At C A 
At CB 

At C B 



162.2+^ 



9.4 
: 23.7 : 

: 25.0 



3.4 



JVAS B 1422+231 



contradictory results 
between methods: 
BACorCAB? 



AtBA ■ 

At A c 
At BC 



1.5: 
7.6: 
8.2: 



-3.j 



1.4 
2.5 
2.0 



Burud (2001 1 
|Hjorth eTaEfl2002f 



Jakobsson et al. ( 2005 ) 
Poindexter et al. (20071 
Ulek& Maoz ( 2003j__ 
Wyrzykowskt e faf||2003| 

Morgan et al. ( 2008a) 

"Schechter et alTl 1997") 
Schechter et al]{TW7} 
BarkanatTrW7T 



Patnaik & Narasimha 



(2001 



SBS 1520+530 


At AB = 125.8 i 


2.1 


At AB 
At AB 


= 130 + 3 
= 130.5 + 2.9 


Burud et al. (2002c ) 


Gaynullma et al. ( 2005b I 


CLASS B 1600+434 


At AB = 47.8 ± 


1.2 


At AB 


= 51+4 


Burud et al. (2000) 


CLASS B 1608+656 


At BA =31.6 + 


1.5 


At BA 


= 31.53 


Fassnacht et al. ( 2002 1 




At BC = 35.7 ± 


1.4 


At BC 


= 36.0+1.5 




At BD = 77.5 ± 


2.2 


At BD 


= 77.0+2 




HE 2149-2745 


unreliable 




At AB 


= 103 ± 12 


Burud et al. (2002a) 



Table 2: Summary of Time Delays for 11 Lensed Systems. 
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